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ABSTRACT 

Context. There have been a long string of efforts to understand the source of the variability observed in microquasars, especially 
concerning the elusive High-Frequency Quasi-Periodic Oscillation. These oscillations are among the fastest phenomena that affect 
matter in the vicinity of stellar black holes and therefore could be used as probes of strong-field general relativity. Nevertheless, no 
model has yet gained wide acceptance. 

Aims. The aim of this article is to investigate the model derived from the occurrence of the Rossby wave instability at the inner edge 
of the accretion disk. In particular, our goal here is to demonstrate the capacity of this instability to modulate the observed flux in 
agreement with the observed results. 

Methods. We use the AMRVAC hydrodynamical code to model the instability in a 3D optically thin disk. The GYOTO ray-tracing 
code is then used to compute the associated light curve. 

Results. We show that the 3D Rossby wave instability is able to modulate the flux well within the observed limits.We highlight that 
2D simulations allow us to obtain the same general characteristics of the light curve as 3D calculations. With the time resolution we 
adopted in this work, three dimensional simulations do not give rise to any new observable features that could be detected by current 
instrumentation or archive data. 

Key words. Accretion, accretion disks - Hydrodynamics - Instabilities - Radiative transfer - Methods: numerical 



1. Introduction 



Black hole binary systems exhibit a wide range of variability 
patterns. Among these, the high frequency quasi-periodic oscil- 
lations (HFQPOs) appear as small narrow peaks in the power 
density spectrum. Their frequencies, of several tens to a few hun- 
dred Hertz, correspond to the rotation frequency in the inner part 
of the disk suiTounding the black hole (see more details in the 
review of ,van der Klis_2006 1 making them an interesting tool to 
study strong-field gravity. Even though they are the fastest phe- 
nomena observed in the vicinity of black holes, HFQPOs are also 
extremely elusive, especially when compared to the lower fre- 
quency quasi-periodic oscillations (LFQPO). For this reason, we 
do not have a rich domain of observations from which to draw 
constraints for the models. This has led to a wide variety of mod- 
els (see e.g. Lai & Tsang|2009" for a very brief review) focusing 
on one or another of the observed characteristics of the HFQPOs. 
As more and more data are obtained, especially with the next 
generation telescopes such as LOFT ( jBozzo et al.|201 1 1 



we will 

be able to gain a better picture of what a model must explain. 
Nevertheless, we already can make a stringent list of constraints 
based on the eight objects with recurrent HFQPOs (see for ex- 
ample the review by Remillard & McClintock 2006 1. Among the 
most important of those constraints we note: 



1- The emitted flux is modulated at a level of a few percent in 
X-ray and this modulation increases with energy. This first 
and fundamental point, often not discussed by models, is the 
aim of this paper. 

2- The frequency of the modulation has a small but significant 
range that must be explained. 

3- The HFQPOs occur sometimes, but not always, in pairs close 
to a 2:3 ratio. Those occurrences, and not just the close ratio, 
need to be elucidated. 

4- Another point often neglected in models for HFQPOs is the 
fact that, sometimes, they co-exist with the more ubiquitous 
LFQPOs. 

The model we are investigating here is based on the exis- 
tence of the Rossby wave instability (RWI) at the inner edge of 



the disk as was discussed in Tagger & Varniere ( 2006 1. Indeed, 



it was shown that the existence of an innermost stable circular 
orbit (ISCO) around a black hole makes the disk prone to the 
RWI because the vorticity profile has a natural extremum. This 
dynamic instability results in the formation of large-scale spiral 
density waves and Rossby waves that can reach high amplitudes. 
Depending on the disk's physical parameters, different modes of 
this instability will be selected, most often with azimuthal mode 
number m = 2 or 3 (Tagg er & Varniere 2006) which gives a 
natural explanation for some of the observed characteristics of 
HFQPOs. Nevertheless, the important step of computing the ac- 
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tual flux modulation one would observe when such an instability 
occurs in a disk was still missing. Here we focus on this point 
and therefore restrict ourselves to the case of the HFQPO alone. 
Indeed, we are trying to explore the capacity of one particular 
HFQPO model to modulate the X-ray flux and having only one 
instability in the disk makes the results more conspicuous. To 
this end we will use the purely hydrodynamic version of the 
Rossby wave instability (RWI, Lovelace et al. 1999), without 
the e xtra destabilizing effect p rovided by a poloidal magnetic 
field ( Tagger & Varniere|2006|l or the ex tra stabilizing effect of 
a toroidal magnetic field ( |Yu & Li|2()09| . 

This is the continuation of the long string of studies on the 
impact of the RWI in disks such as was done in the context of 
protoplanetary disks ( Regaly et al. |20I2 1 and the galactic cen- 
tre ( Falanga et al.|2007 1. Whereas the first case is highly differ- 
ent from ours, the second is more comparable and we use sim- 
ilar tools, but Falanga et al. (2007i restricted themselves to the 



2D case and to the special context of the galactic centre source 
Sgr A*. 

Our main interest here will thus be to answer the key ques- 
tion: is the RWI able to modulate the flux of a microquasar at 
a level coherent with observations? In a more practical manner, 
we will also address the question of the necessity of full 3D sim- 
ulations of the RWI, and to what extent it is sufficient to restrict 
oneself to 2D simulations. Those 2D simulations are more acces- 
sible in terms of computing and storage resources, therefore al- 
lowing us to do much larger parameter studies and with a higher 
time resolution that could be compared with future LOFT ob- 
servations. These questions will be addressed by computing the 
light curve of a microquasar subject to the RWI, as observed by a 
distant observer. This computation will be performed by means 
of a general relativistic ray-tracing code ( [Vincent et al.|[20TT] l. 
The main general relativistic effect that affects the RWI being 
the existence of an ISCO, we use a pseudo-Newtonian potential 
in order to model this in the hydrodynamical simulations, instead 
of considering the full general relativistic case. This simplifying 
assumption is sufficient in order to derive a proof of principle of 
the ability of the RWI to modulate the flux of microquasars. 

Sect. |2] briefly presents the RWI and then develops the hy- 
drodynamical simulations that allow us to model the unstable 
accretion structure. The following section is focused on the ray- 
tracing procedure used to compute observable quantities. Sect.|4] 
presents the light curves that were obtained and Sect. |5] gives 
conclusions and prospects. 



2. Hydrodynamical simulation of the RWI 

The RWI has been discussed in multiple astrophysical contexts, 
from the galactic centre to protoplanetary disks. It can be seen as 
the form the Kelvin-Helmholtz instability takes in differentially 
rotating disks, and has a similar instability criterion. For two- 
dimensional (vertically integrated) barotropic disks, the RWI can 
be triggered provided an extremum exists in a quantity X which 
is the inverse vortensitjQ 



SQ p 

2^ET'' 



(1) 



where p is the pressure, S is the surface density, Q is the rotation 
frequency, — 2f2/rd/dr (r^Q) is the squared epicyclic fre- 
quency and J is the adiabatic index. An extremum of vortensity 



can thus typically arise from an extremum of the epicyclic fre- 
quency, as is the case close to the ISCO, or from an extremum 
in the density profile which gives rise to the RWI in the disk. 
Once it is triggered, it leads to large scale spiral density waves 
and Rossby vortices. The dominant mode (i.e. the number of 



vortices) is dependent on the disk conditions (Tagger & Varniere 
[20061 ). 

2. 1 . The RWI as a model for the HFQPO 

Due to general relativistic effects, the epicyclic frequency profile 
will show a maximum in the vicinity of the black hole's ISCO, 
creating a vortensity extremum . [Tagger & Vamierel ( [2006| l 
showed that it would lead to the growth of the RWI at the inner 
edge of black hole binaries when the disk nears its ISCO. From 
these facts, it seems natural that the RWI model may be capable 
of accounting for points 2- to 4- in the Introduction above (that 
will not be investigated in detail here). 

2- The frequency range: the exact location of the resulting 
vortensity extremum will depend on the density profile of 
the disk on top of the epicyclic frequency. Thus the location 
of the launching of the instability can vary depending on the 
disk properties. As the Rossby vortices will orbit around the 
central black hole with the disk's rotation frequency at the ra- 
dius of launching, this variation of the disk properties gives 
rise to a range of different possible frequencies of the signal. 

3- Pairs of frequencies: different modes can dominate the sig- 
nal depending on the physical state of the disk. This was 
ad ressed in [Tagger & Varniere| ( |2006| and more recently 



Varniere et al. ( [2012b i. The 2:3 frequency pairs can be 



naturally explained by a superposition of azimuthal modes 
m — 2 and m — 3. 
4- Coexistence of HF- and LFQPO: we have tackled this point 
by showing the ability of the RWI to co-exist with an- 
other instability that we proposed to be at the origin of 
the LFQPO ( Varniere et al.[[20I2a| l, the Accretion Ejection 
Instability (AEI). A disk giving rise to both AEI and RWI 
will thus naturally show both types of QPOs. 

It seems therefore natural to address in more detail if the 
RWI is strong enough to give rise to a modulation that could 
explain the observations. 



2.2. Numerical setup 

To perform the hydrodynamical simulation of the RWI we use 
the Message Passing Interface-Adaptive Mesh Refinement 
Versatil e Advectio n Cod e (MPI-AMRVAC) developed 
by Keppens et al. ( 201 l| l. The numerical scheme is the 
le for all refinement levels, namely the Total Variation 



all 



Diminishing Lax-Friedrich scheme (see Toth & Odstrcil[[T996| l 
with a third order accurate Koren limiter ( [Koren[[1993[ ) on the 
primitive variables. 

We consider a geometrically thin disk and neglect self- 
gravity. The Euler equations in cylindrical coordinates (r, tp, z) 
read: 



d,p + V ■ (vp) = , 
d,(pv) + V ■ (vpv) H- Vp = -pV^G : 



(2) 
(3) 



Vortensity is the ratio of vertical vorticity to surface density. 



where p is the mass density of the fluid, v its velocity, and p 
its pressure. We consider a barotropic flow, i.e. the entropy is 
constant in the entire system. The pressure is then p = Sp'^, 
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with the adiabatic index j - 5/3 and the constant 5 related 
to entropy. The sound speed is given by c\ - ypip - Syp^^^ 
and the temperature by T oc pip oc Sp'^^\ All the distances 
are normalized by the ISCO radius ry. The same applies for 
temperatures normalized at ISCO by: To = 10^ K, while all 
times are in units of the ISCO period fisco- 

The general relativistic effects have to be taken into account 
to study the innermost region of the black hole disk. In order 
to provide a proof of principle of the ability of the 3D RWl to 
modulate the flux of microquasars, we use the minimum model 
necessary to mimic some of these general relativistic effects by 
using a pseudo-Newtonian potential Oc - For our simulations, we 
use the common Paczynsky & Wiita ( 1980[ l gravitational poten- 
tial: 



1.2 



1.4 



1.6 



2.0 



r/ro 



Of; = -- 



where Kg - is the gravitational radius, G the gravitational 
constant, M» the black hole mass, c the speed of light in vac- 
uum. With this gravitational potential, one can define the inner 
limit below which there is no stable circular orbit, the ISCO, lo- 
cated at rg = 6 Tg. Let us stress the fact that this choice of poten- 
tial allows us to mimic only some aspects of the Schwarzschild 
metric. The effect of the black hole spin is not taken into account. 
It is also important to keep in mind that this pseudo-Newtonian 
potential, although widely used, is far from giving an exact de- 
scription of general relativity. However, the only characteristic 
of the Schwarzschild metric that is of primordial importance as 
far as RWl is concerned is the existence of an ISCO, and this is 
correctly modeled by Eq. |4] 

Let us mention that we have also developed 2D RWl sim- 
ulations with alternative pseudo-Newtonian potentials that take 
into account the rotation of the black hole. These alternative po- 
tentials were taken from I Art emova et"aL] ( |1996| l. The RWl also ifi z) - Pm( 1 - 
develops well in these potentials, and future works will be ded- \ 
icated to studying the effect of the black hole's rotation on the 
observables. 



2 r„ 



(4) 



Fig. 1. Inverse vortensity profile, normalized to its extremum 
value. The instability is triggered at the location of the ex- 
tremum, i.e. at r » l-4ro. 



rms of the modulation depends on this choice of slope. However, 
we are here only interested in proving the existence of the flux 
modulation, not comparing a precise value of modulation to ob- 
servations. If future instruments give better data on the rise of 
the HFQPO we might be able to use these to constrain the disk 
model that could reproduce the observed data. 
The initial density profile then reads: 



Pm = p(r, (p,z^O) 



Po 
2 



the hydrostatic equilibrium gives the vertical profile of the 
density 



r-1 



ySp 



M 



GM, 



r-2rg Vr2 +z^-2r. 



(6) 



2.3. Disk setup and boundary conditions 

The grid is cylindrical with the unit length defined as the ISCO 
radius tq. The radial coordinate r is in the range [0.8, 6], the full 
azimuthal coordinate spans [0, 27r] and the vertical coordinate z 
lies in the range [0,0.2] for the 3D simulations. For the fluid 
simulations, we considered only the upper part of the disk as 
the mid-plane is a symmetry plane for the RWl (Meheut et al. 
[20T0l 



2012 



. For the ray-tracing computations, the full disk is 
)y symmetrizing the upper part. We used a fixed and 

(384,72) in 2D 



considered 

homogeneous grid with a resolution rir x 



and X x n- - (384, 72, 72) in 3D, which is slightly higher 
than the resolution used in Meheut et al.| ( |2010| . 

As we describe gravitation by the pseudo-Newtonian poten- 
tial in Eq. |4] the epicyclic frequency will show a maximum at a 
radius slightly higher than the ISCO radius tq. This epicyclic 
frequency maximum gives rise to an extremum of the initial 
vortensity profile, which is shown in Fig. [T] The existence of 
this vortensity extremum does not depend on the choice of ini- 
tial density profile even though it does influence its precise lo- 
cation in the disk. Considering the density, we choose a typical 
power-law profile with a slope of -3/2. This choice, while rea- 
sonable, allows us to obtain a steeper vortensity extremum and 
hence a stronger instability. Note that the precise value of the 



where rg - 1.3 gives the position of the density maximum (its 
value was choosen to fit the simulations of Tagger & Varniere 
|2006 |, rg - ro/6 is the gravitational radius, cr - 0.05 ro gives the 
width of the plunging region, a - -3/2 is the density slope, and 
Po is the minimum density of the simulation. These parameters, 
and mainly the choice of density slope, can modify the structure 
and characteristics of the RWL In our case, the existence of 
an extremum in k limits the impact of the density power-law 
index and any reasonable choice will exhibit the instability. 
The detailed effect of the initial density profile on the growth 
rate of the instability and its saturation level will be studied in 
a forthcoming paper (Meheut, Lovelace, Lai, 2013) but will not 
influence its ability to modulate the flux. 

We also considered a 2D disk with the surface density de- 
fined as S oc r"^^^ in the outer region. The initial conditions 
of the 2D and 3D simulations are then highly different. In both 
cases, the azimuthal velocity is determined by force balance: 



\(r-2r,)2 



+ o:Sypl^^ +Syplj^ 



2crcosh2 



(7) 



The density and velocity profiles given in Eqs. 6]|7 include 
high amplitude gradients and do not give exact numerical equi- 
librium. For this reason, the initial conditions of our simulations 
are numerically computed from this pseudo-equilibrium. This 
means that a first simulation is done without any perturbations 
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5 







Fig. 2. Midplane density map of the 3D simulation when the 
m=l mode dominates. The axes are in units of the ISCO radius. 
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and is run until the disk has reached a permanent stage which is 
chosen as the initial conditions. This disk is then perturbed with 
small amplitude (~ 10"^ v^) random perturbations on the radial 
velocity which act as seeds for the instability. 

The inner boundary condition is a no-inflow condition and 
the outer boundary condition is a null radial velocity condition. 
Meheut et al. (2012 1 have shown that the boundary conditions 
do not change significantly the growth rate of the RWI due to its 
confinement in the vortensity bump. Moreover the inner edge of 
the simulation being inside the ISCO, any reflected wave would 
not reach the region of interest for the instability. Therefore, the 
boundary conditions are not determinant for these simulations. 
For the 3D simulations, symmetric boundary conditions are im- 
plemented at the mid-plane, and we use a null vertical velocity 
at the grid upper boundary limit, situated outside of the disk. 

2.4. Time evolution of the RWI 

Fig. |2] shows the density map of the z = plane of a 3D RWI 
simulation, when the instability is completely developed. This 
section presents the time evolution of the instability (at 2D and 
3D) from its launch to its complete development. 

At first, the disk is assumed to be 2D and the vertical struc- 
ture of the disk is neglected. The growth of the instability is 
shown on Fig. [3]where the time evolution of the perturbation of 
density is plotted on a logarithmic scale. This allows the identi- 
fication of the linear phase of the instability when the amplitude 
of the perturbations grows exponentially (thus linearly in loga- 
rithm), as well as the saturation which is due to non-linearities. 
In the 2D case, the initial conditions give rise to a rapid growth 
of the instability on a time scale of a few periods at ISCO as can 
be seen in Fig.|3] During the linear phase, the perturbations can 
be separated in the different azimuthal modes: 



X 



^ X,„ sin(mwf ■ 

meN 



m4)) expy„,f 



(8) 



where X - 1. for 2D simulations and X - p for 3D. The 
quantities y„, are the growth rates of each mode, u the charac- 
teristic frequency of the instability depending on the position 
of the extremum of £,, m the azimuthal mode number, and X,n 
the amplitude of mode m. Therefore, X,„=o corresponds to the 
axisymmetric part of the density and the frequency of the mode 
m is the multiple mu of the frequency of the fundamental mode. 
Fig. [3] also shows the time evolution of the amplitude of the 



Fig. 3. Amplitude of the surface density perturbation in logarith- 
mic scale (solid line) for the 2D simulation. The amplitudes of 
the strongest modes are also plotted in dotted and dashed lines. 
The vertical dashed lines show the changes of dominant modes, 
indicated in italic at the top of the figure. The fundamental mode 
m - \ eventually dominates. 




t/t,s, 



Fig. 4. Same as Fig. [3] for 3D simulations. However, mind the 
different time scale. 



strongest modes. During the linear phase, the dominant mode is 
m - 3, and later on the disk is dominated by the fundamental 
mode m - I with important contributions from the modes m - 2 
and 3. This evolution of the oscillation modes depends on the 
disk's astrophysical properties: different modes dominate for 
different initial disk conditions. 

This fact is illustrated in the 3D simulations. Indeed, the ini- 
tial conditions differ largely between the 2D and 3D simulations, 
with a different surface density radial profile and absolute value. 
The characteristic velocity of the waves is then different. This 
modifies the timescale for the growth of the instability as can 
be se en in Fig. |4l T his difference is not due to the dimension- 
aUty ( [Meheut et aL][20TH pnl[20T2| . Nevertheless, the RWI is 
triggered and saturation of the 3D instability is reached after a 
few tens of periods at ISCO. The dominant modes are also differ- 
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ent during the linear phase, but non-linearities still tend to favour 
the lowest azimuthal mode. 

Let us note that, as the epicyclic frequency profile is not 
evolving with time, the location of the extremum of vortensity 
will stay the same. Moreover the vortices grow in a density ex- 
tremum where they cannot migrate due to the two density slopes 
of opposite sign. As a consequence, the Rossby vortices will not 
migrate during the simulation (for more details see Meheut et al. 

From the hydrodynamical simulations we obtain the density 
and velocity at every point of our grid with a sufficient time sam- 
pling (of 40 frames in 2D and 6 frames in 3D per orbit at the 
ISCO) to follow the RWI in the disk from its rise to its saturation. 
In the next step those will be used as input for the ray-tracing 
code in order to trace the impact of the RWI on the observed 
emission. 

3. Ray-tracing of the RWI with GYOTO 

3.1. Ray-tracing the accretion disk 

In order to obtain images, and especially lightcurves, that could 
be compared with observations we used the general relativistic 
ray-tracing code GYOTO (Vincent et al. 2011). Null geodesies 
are computed backward in time in the Schwarzschild metric 
from a distant observer to the emitting disk. The coordinates 
used by GYOTO are spherical-like and denoted (r, 0, (f), where 
r is used to diff'erentiate the GYOTO spherical radius from the 
fluid simulations cylindrical radius r. Once a backward inte- 
grated geodesic hits the disk, the density and 3-velocity at the 
point of emission is linearly interpolated at the time of emission 
from the results of the MPI-AMRVAC computations. 

3.1 .1 . Radiative processes at 2D: blackbody 

The 2D disk is assumed to be optically thick, so the integration 
is stopped at the first encounter of the disk. The emission is sup- 
posed to follow the Planck law; the only parameter needed is 
thus the temperature which can be easily derived from the den- 
sity according to the following computations. 
The gas being assumed ideal: 



P 



P 



kT 



(9) 



where jj is the mean molecular weight, ;«„ is the atomic mass 
constant and k is the Boltzmann constant. Assuming that the 
plasma is made of pure hydrogen, = 1. Using the relation 
R - NAk between the ideal gas contant R, the Avogadro number 
Na and the Boltzmann constant together with the expression of 
the sound speed this yields: 



Na niu 

R 



7 



(10) 



This allows the computation of the emitted specific intensity 
at the surface of a 2D disk: 

ly^BJT) (11) 

where By is the Planck function. 

3.1.2. Radiative processes at 3D: Bremsstralilung 

For the three dimensional computations, the only radiative pro- 
cess considered is Bremsstrahlung. As the disk is purely hydro- 
dynamic, there is no synchrotron emission, and Compton scat- 
tering is neglected as we are only interested here in a proof of 



principle, not in a detailed study of emission processes. The 
Bremsstrahlung emission is assumed to be thermal, so that the 
emission coefficient jy and absorption coefficient ay are related 
via Kirchhoff's law: 



jy - ay By 



(12) 



The emission coeflScient for thermal Bremsstrahlung is given 
by ( [Rybicki & Lightman|1979l ): 



1 l^ne^ I In 



Jy 



An 3mpC^ \ 3km, 



1/2 



,-1/2 rp 



exp 



hv 
"kf 



(13) 



where e is the electron charge, me is the electron mass and h is 
the Planck constant. Here, we assume that the disk is made of 
pure hydrogen, and that the emission is isotropic in the emitter's 
frame (hence the 1 /An initial factor). Moreover, the Gaunt factor 
is neglected as most of the radiation arises from locations in the 
disk where hv x kT. 

Once the emission coefficient is computed, the absorption 
coefficient is derived by using Eq. 12 



3.1 .3. Dependency on temperature at 2D and 3D 

Let us investigate the dependency on temperature of the 3D 
emission process (Bremsstrahlung) as compared to the 2D case 
(blackbody). For temperatures around lO^K, from where most 
of the flux arises, the 3D emission coefficient follows jy'' oc 
T^-^exp(-hv/kT) whereas the 2D specific intensity follows 
7*^ oc exp(-hv/kT). The dependency on temperature is thus 
much more important in 3D, and the emission arises only from 
the hottest parts of the disk, whereas the emission is more spread 
out in the 2D case. 

All these computations allow us to determine the specific in- 
tensity in the emitter's frame, ly^^. In order to compute the spe- 
cific intensity in the observer's frame ly^^^^, the frame invariant 
quantity ly/v^ is used: ly^^^ = '•'obs/^em ^v„,- The quantity Vem 
can be related to the emitter's 4-velocity. The 3-velocity com- 
puted by the fluid simulations is used to determine the emitter's 
4-velocity, and the observed specific intensity is thus at hand. 



3.2. Computing tine disk image and tine iigint curve 

Before producing an image of the disk, one needs to compute 
the relativistic time delay for each point of the disk in order to 
take into account the multiple time steps of the fluid simulations 
that will correspond to the same observed time. Fig.|5]shows the 
emission date of photons for two extreme inclinations (5° and 
85°). Fig.|5]shows that the emission points are spread differently 
along the disk depending on the inclination. The distribution 
is homogeneous and isotropic for low inclinations, but is very 
inhomogeneous and anisotropic for high inclinations. This 
anisotropy is due to the projection effect due to the mapping of 
the observer's screen onto a very inclined disk, resulting in the 
emission points being aligned preferentially along the direction 
perpendicular to the line of sight. The inhomogeneity is due to 
the tensing effect of the black hole that concentrates emission 
points on the side of the black hole opposite to the observer. 
Fig.|5]also clearly shows that the higher the inclination, the more 
time slices of data will be required. The respective difference 
between the maximal and minimal emission dates on the pri- 
mary image are approximately 0. 1 fisco (left panel) and 0.9 fisco 
(right panel). However, due to the strong beaming effect at high 
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observer 
> 



-0.45T,sco -0.45T|sco 



Fig. 5. Emission dates of photons at the footpoint of null geodesies connected to the observer's screen, plotted with the same color 
scale for an inclination of 5° (left) and 85° (right). The grid of emission positions is projected onto the plane of the accretion disk. 
The color bar is expressed in units of the ISCO orbital time fisco- The lensing effect of the black hole is clear for an inclination of 
85° on the right panel, leading to a higher concentration of emission points on the side of the black hole opposite to the observer. 



inclination, the bulk of the total specific intensity in one image 
comes from a small part of the disk. Thus, the emission dates 
of the photons reaching this small part of the disk are close to 
each other, and the final effect of time delay is only marginal on 
the light curve. We have checked that the difference between the 
exact light curve and a light curve computed without taking into 
account the time delay is at the level of one to a few percent only. 

Fig.|6]shows the image (i.e. the map of specific intensity) of 
the accretion disk at an inclinatiorj^of 85°, approximately at the 
time when the fundamental mode of the RWI dominates. Each 
pixel of the image is obtained by interpolating between the set of 
simulated data at the time of emission of the photon, as stressed 
in Sect. [3] Here, around 40 different data time slices are used for 
computing the image. The instability has started to grow and one 
can identify the spiral density waves of the RWI. On top of this 
modulation due to the instability, there is a clear beaming effect: 
matter moving towards the observer is brighter (here, on the left 
side of the image). The secondary image of the disk is visible as 
a semicircle at the centre of the primary image. 

Fig.|7]depicts the image of a 3D disk subject to the RWI, ap- 
proximately at the time when the fundamental mode of the RWI 
dominates. Here, around 6 different data time slices are used for 
computing the image: the time sampling has been reduced as 3D 
data are much heavier than 2D data. We checked at 2D that this 
reduced sampling still allows us to retrieve similar results. As 
compared to Fig. |6j the 3D disk's emission is much more con- 
centrated on the inner parts of the disk. This is due to the much 
stronger dependency of the 3D emission on temperature as com- 
pared to the 2D case, as explained in Sect. 3.1.3 As the hottest 




parts of the disk are concentrated close to the radius of launching 
of the instabihty (see Fig.|2]|, only these inner regions are seen. 



Fig. 6. Image of a 2D optically thick disk emitting blackbody 
radiation with an inclination of 85° and in the presence of the 
RWI . The first order image (the distorted complete ring) clearly 
shows the spiral shape of the emitting region. The beaming effect 
makes the emission brighter on the approaching side of the disk 
(here, on the left of the image). The second order image is the 
portion of ring at the bottom of the image. It is due to photon 
swirling around the black hole before reaching the observer. The 
third order image is the thin ring of illuminated pixels at the 
center of the image, due to photons orbiting around the black 
hole very close to the event horizon before reaching the observer. 
This ring is the so-called black hole silhouette, it is truncated 
here due to optical thickness. 



4. Modulation of the light curve by the RWI 

This section presents and analyzes the light curves obtained 
by ray-tracing the RWI hydrodynamical simulations. These are 
computed by summing the specific intensities over all solid an- 



^ We call here inclination the angle between the z axis and the line of 
sight. A null inclination thus corresponds to face-on view. 



gles on the observer's sky. This boils down to summing the disk 
images over all pixels, for all observation times. 

4.1. 2D case: high time resolution 

We first analyze the main characteristics of the light curve 
modulation in the 2D case where we can have a much better 
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Fig. 7. Image of a 3D disk emitting Bremsstrahlung radiation 
with an inclination of 85° and in the presence of the RWI. The 
overall emission is much more concentrated in the vicinity of the 
Rossby vortex as compared to Fig.|6] 

time resolution between the hydrodynamical snapshot^ In 
order to catch the details of the light curve evolution, 40 frames 
of fluid simulations are computed during each period of rotation 
at ISCO. The ray-tracing code then interpolates between these 
fluid simulation results in order to compute the specific intensity 
map at dififerent times of observation. We compared the results 
obtained with the one from a higher (80 frames) and smaller 
(20 frames and 6 frames) time resolution. While the overall 
behaviour was similar, 40 frames allow us to get a more detailed 
light curve, similar to that with 80 frames. We therefore decided 
to do all the runs in 2D at 40 frames per period of rotation at 
ISCO. We computed the light curve from the moment the RWI 
started to grow in the disk until it reaches its non-linear states as 
shown by the amplitude evolution in Fig. |3] 

At zero inclination, the light curve will slowly drift towards 
smaller values due to the black hole's accretion during the sim- 
ulation. In order to determine the flux fluctuation that is only 
due to the RWI, it is important to subtract this continuous drift 
from the light curve. This is done by simply subtracting the light 
curve obtained at 1° of inclination from every other light curve. 
As the GYOTO code uses Boyer-Lindquist spherical-like coor- 
dinates, the z-axis is singular, and it is not possible to derive a 
light curve at exactly zero inclination. However, the light curve 
is only marginally impacted by beaming at 1° as compared to ex- 
actly 0° . The error introduced by subtracting the light curve at 1 ° 
of inclination is a small underestimation of the amplitude of the 
modulation that does not impact our result here. Indeed, we are 
only looking for a proof of principle that the RWI could modu- 
late the flux at a level compatible with the observed HFQPO. 

Fig.[8]shows the light curves obtained for different values of 
the inclination parameter at two different values of the energy 
of the observed photon: 1 and 2 keV. Initially when the insta- 
bility has a very low amplitude, the modulation of the flux is 
very low. After a few ISCO rotations, the instability is growing 
exponentially and modulates the flux at a detectable level. The 
period Tmod of the modulation equals the period of rotation of 



^ This is due to the smaller size of the 2D data files compared to the 
3D. 



the Rossby vortices, i.e. the period of rotation at r s; 1.4 ro (see 
Fig.[T}, that is to say T^^a ~ 1-6 fisco- 

Due to the general relativistic beaming effect, when the 
Rossby vortex is on the approaching side of the disk, the result- 
ing flux is boosted at high inclination, the opposite being true 
for the receding side of the disk. This beaming effect is also vis- 
ible in the light curve substructures that can be observed at high 
inclination. For instance, an m - 2 mode will lead to two sharp 
peaks in the light curve at high inclination, related to the passage 
of the Rossby vortices at the approaching side of the disk. If the 
inclination is low, this effect is much fainter, leading to smaller 
substructures in the light curve. The oscillation frequency of the 
light curve evolves during the simulation with a high frequency 
after a few rotations. After 5 rotations, a mixture of modes be- 
tween one frequency and twice this frequency, can be identified. 
Whereas at the end of the simulation the mode with the lowest 
frequency dominates. This evolution corresponds to the evolu- 
tion of modes seen in Fig. |3] the dominant mode has initially 
a frequency 3a>, then 2a) and eventually the fundamental mode 
dominates. 

The comparison between the two values of energy of the ob- 
served photon shows that similar shapes are obtained for the two 
cases, but a higher modulation amplitude is reached for higher 
energy as expected due to the Planck law. Indeed, the ISCO tem- 
perature being 10^ K, the maximum frequency of the Planck law 
is closer to 2 keV that to 1 keV. 

Fig. |9] shows the light curve rms computed with a sliding 
window with a width of two orbital periods. The amplitude of 
the modulation is growing linearly with a change of slope at 
around 5 fisco, corresponding to the saturation of the instabil- 
ity (see Fig.[3]l. The maximum level of modulation is of around 
4% at 1 keV and 8% at 2 keV. This is somewhat higher than 
most observations of HFQPOs. However, the maximum's exact 
value is influenced by the initial X. profile (see Eq.[T and Fig.[T]i. 
Here we are only interested in seeing if a reasonab e setup can 
modulate the flux to a level similar to the one observed. We keep 
the detailed study of the growth rate and saturation level for a 
forthcoming paper (Meheut, Lovelace, Lai, 2013) 

4.2. Full 3D case 

We then turned to 3D simulation to see if there was any differ- 
ence in the observables. There, we could not get a similar time 
resolution as in 2D because of the hydrodynamical data size. We 
therefore used the 40 frames per orbit resolution for only one 
period for confirmation while we made a longer term lightcurve 
at the much lower time resolution of 6 frames per orbit at the 
ISCO. As we will see below, this lower resolution is still able 
to catch the broad aspects of the light curve, although it will not 
allow us to obtain the finest details. 

Fig.[TO]depicts the light curve and its rms for a 3D disk sub- 
ject to the RWI, for photons of observed energy equal to 2 keV, 
and for three different values of inclination. The first important 
result is that the RWI is capable of modulating the light curve, 
similar to what had already been shown in the 2D case, with 
a modulation of a few percent. This is a strong argument that 
makes the RWI a reliable model of microquasars HFQPOs. 

The domination of the different modes is also easily seen in 
Fig. 10 with the mode m = 2 dominating at the beginning of 
the simulation, and eventually the mode m = 1 . Here, the mode 
m = 1 dominates after around 15 ISCO periods, whereas it dom- 
inates after around 30 periods in Fig. |4] This difference is due 
to the fact that the ray-tracing simulations are initiated when the 
instability is already strong enough to give a non negligible rms. 
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Fig. 9. Dispersion of the light curve points in Fig. |8jat 1 keV (left) and 2 keV (right). Each point is obtained by computing the 
dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds), or 85° (red 
crosses). 



As a consequence, the 15 first ISCO periods of the 3D hydrody- 
namics computations where not used as their rms is extremely 
low. 

The right panel of Fig. [TO] shows the same characteristics as 
in the equivalent 2D Fig. [9] The rms shows a linear profile with a 
clear change of slope around f ^ 7 fisco- This is linked to mode 
m - I starting to dominate over mode m = 2. 

This being given, there are two main differences between the 
2D and 3D results: 

1- The 3D modulation grows more slowly than its 2D counter- 
part. 

This is comes from the differences in initial disk condi- 
tions, i.e. mainly to the choice of the density profile: differ- 
ent density profiles give different shapes for the extremum 
of X which in turns give different growths of the instabil- 
ity. Figs. [3] and |4] already showed that the 3D growth rate is 
slower than its 2D counterpart. This translates directly to the 
light curve (moreover, as stated above, the 3D light curve 



initial point is 15 ISCO periods after the launch of the insta- 
bility). 

Let us stress that the index of the density profile is not con- 
strained by observations. As the difference of growth time 
depends on the choice of the power-law index, it cannot not 
be taken as an intrinsic difference between 2D and 3D RWI. 
The light curve exhibits a slight asymmetry in the 3D case 
(the oscillation is not exactly centered on the initial flux 
value). 

This could be explained by the strong dependency on tem- 
perature of the 3D emission, as explained in Sect. 3.1.3 



During the simulation, the temperature of the Rossby vor- 
tex increases by a few %, due to accretion. Its emission thus 
becomes greater and greater. This translates to a slight drift 
of the light curve maxima towards higher values. 
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Fig. 10. Left: Light curves of a microquasar subject to the 3D RWI at 2 keV at an incUnation of 5° (dash-dotted black), 45° (dashed 
blue) or 85° (solid red). Right: Dispersion of the light curve points of the left panel. Each point is obtained by computing the 
dispersion of the light curve points over two orbital periods. The inclination is 5° (black triangles), 45° (blue diamonds) or 85° (red 
crosses). For both panels, mind the different time scale as compared to Figs. [8] and |9] 



4.3. Are 3D simulations necessary to compare with 
observations? 

Comparing Figs. |8]to [TO] it appears first that the light curve is 
modulated both in 2D and 3D, which is the main result of this 
paper. The dependency of the light curve as a function of the in- 
clination parameter and frequency of the radiation is also very 
similar. As these general features of the 2D and 3D computa- 
tions are alike, an interesting conclusion is that 2D results are 
sufficient in order to analyze the general observable characteris- 
tics of the RWI, at least until we get a better constraint on the 
profiles in the disk. 

This is particularly interesting when considering the differ- 
ence in terms of computing time and memory resources between 
a 2D and 3D simulation. The hydrodynamical data at a given 
time is around 100 times heavier at 3D than 2D (typically re- 
spectively 50 Mo and 0.5 Mo). The time needed to compute an 
image from 3D data is typically 15 times what is needed for a 
2D computation. 

However, let us stress that the present work does not allow us 
to compare the detailed time evolution of the 3D RWI with the 
2D case, due to our choice of time sampling in the ray-tracing 
computations (this choice being dictated by the computing time 
and memory resources needed for the 3D simulations). As the 



3D RWI displays specific features in the z direction (Meheut 
et al.| 2010 1, it may be that specific observable characteristics 
could be obtained only by resorting to high time resolution 3D 
simulations. Nevertheless, these features would be small correc- 
tions to the general trend that stays close to the 2D results, and 
would not be within reach of current instruments. Indeed the 2D 
time sampling of 40 frames per ISCO period implies a time res- 
olution of around 10 s for the observed light curve. This is far 
beyond current instrumental capabilities. 

Moreover, let us stress that the radiative processes and radia- 
tive transfer are treated in a very simplified way in this study for 
the 3D case. The 2D simulations are thus sufficient only when 
one is not interested in studying precisely the radiative proper- 
ties of the disk. 



5. Conclusion 

The RWI has been previously proposed as a model for HFQPOs 
and we have now demonstrated its ability to modulate the flux 
coming from the disk. Using 2D and 3D hydrodynamical simu- 
lations we have also been able to study how the amplitude of this 
modulation evolves as a function of the source inclination and of 
the radiation frequency. The 2D simulations have been shown to 
be sufficient in order to recover the broad characteristics of the 
light curve. 

By using hydrodynamical simulations we were able to fo- 
cus only on the RWI and its effects on the light curve. If this 
can be considered for the case where HFQPOs occur alone in 
a softer state, as is the case for the 67Hz modulation of GRS 
1915+105 for example (Mor gan et al.|1997| l, HFQPOs are more 
commonly observed in the steep power law or hard intermediate 
state simultaneously with a LFQPO. In Varniere et al. ( 2012a| l 
we have demonstrated, in 2D, the ability of the RWI to co-exist 
with another instability that could give rise to the LFQPO. Our 
future work will be devoted to that particular case and in partic- 
ular how it influences the observables. Indeed, this state is much 
more frequent during microquasar outbursts than the softer state 
we studied here. 

Acknowledgements. This work has been financially supported by the French 
GdR PCHE and Campus Spatial Paris Diderot. Some of the simulations were 
performed using HFC resources from GENCI-CINES (Grant 2012046810). 



References 

Artemova, I. V., Bjoernsson, G., & Novikov, I. D. 1996, ApJ, 461, 565 

Bozzo, E., den Herder, J. W., Feroci, M., & Stella, L. 2011, in Extreme and 

Vaiiable High Energy Sky (Extremesky 201 1) 
Falanga, M., Melia, F, Tagger, M., Goldwurm, A., & Belanger, G. 2007, ApJ, 

662, L15 

Keppens, R., Meliani, Z., van Marie, A., et al. 2011, Journal of Computational 

Physics, doi:10.1016/j.jcp.2011.01.020, 
Koren, B. 1993, A robust upwind discretization method for advection, diffusion 

and source terms. Vol. 45, Notes on numerical fluid mechanics, ed. C. B. 

Vreugdenhil & B. Koren, 117 
Lai, D. & Tsang, D. 2009, MNRAS, 393, 979 
Lin, M.-K. 2012, ArXiv e-prints 

Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. R 1999, ApJ, 513, 805 



9 



F. H. Vincent et al.: Rossby wave instability: toward a HFQPO model 

Meheut, H., Casse, R, Varniere, P., & Tagger, M. 2010, A&A, 516, A31+ 

Meheut, H., Keppens, R., Casse, E, & Benz, W. 2012, A&A 

Meheut, H., Yu, C, & Lai, D. 2012, MNRAS, 2748 

Morgan, E. H., Remillard, R. A., & Greiner, J. 1997, ApJ, 482, 993 

Paczynsky, B. & Wiita, P J. 1980, A&A, 88, 23 

Regaly, Z., Juhasz, A., Sandor, Z., & Dullemond, C. P 2012, MNRAS, 419, 
1701 

Remillard, R. A. & McClintock, J. E. 2006, ARA&A, 44, 49 

Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysics 

Tagger, M. & Varniere, P 2006, ApJ, 652, 1457 

Toth, G. & Odstrcil, D. 1996, J. Comput. Phys., 128, 82 

van der Klis, M. 2006, Rapid X-ray Variability (Compact stellar X-ray sources), 
39-112 

Varniere, P., Tagger, M., & Rodriguez, J. 2012a 

Varniere, P, Tagger, M., Vincent, E H., & Meheut, H. 2012b, in SE2A-2012: 
Proceedings of the Annual meeting of the Erench Society of Astronomy and 
Astrophysics 

Vincent, E. H., Paumard, T, Gourgoulhon, E., & Perrin, G. 2011, Classical and 

Quantum Gravity, 28, 22501 1 
Yu, C. & Li, H. 2009, ApJ, 702, 75 



10 



